#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <omp.h>
#include <math.h>

#define max(x, y) (x > y ? x : y)
#define abs(x) (x >= 0 ? x : -x)


int main(int argc, char** argv) {
  
  
    FILE *ifp; /* Variavel que contem o ponteiro para arquivo de entrada */
    
    /* Verifica se foi possivel abrir o arquivo */
    if((ifp = fopen(argv[1], "r")) == NULL)
    {
        printf("Could not open the file %s\n", argv[1]);
    }

    int order = 0; /* Variavel que indica a ordem da matriz */
    int row_test = 0; /* Variavel que indica a linha a ser avaliada */
    int i; /* Variavel contadora */
    int j; /* Variavel contadora */
    int ite_max = 0; /* Variavel que indica o numero maximo de iteracoes */

    float error = 0.0f; /* Variavel que indica erro para a condicao de parada */

     /* Recebe valores de definicao da matriz */
    fscanf(ifp, "%d\n%d\n%f\n%d", &order, &row_test, &error, &ite_max);
    
    float* Matrix_A = (float*) malloc(order*order*sizeof(float)); /* Declaracao da matriz A usada na CPU */
    float* Vector_x = (float*) malloc(order*sizeof(float)); /* Declaracao do vetor x usado na CPU */
    float* Vector_b = (float*) malloc(order*sizeof(float)); /* Declaracao do vetor b usado na CPU */
    float* Result_CPU = (float*) malloc(order*sizeof(float)); /* Declaracao da vetor resultado na CPU */
//     int*   Iter_CPU = (int*) malloc(order*sizeof(int)); /* Declaracao usado para saber as iteracoes */
    
    float* Matrix_A_e = malloc(order*order*sizeof(float));
    float* Vector_b_e = malloc(order*sizeof(float));
    
    /* Recebe os valores do arquivo para a matriz A */
    for(i = 0; i < order; i++)
        for(j = 0; j < order; j++)
            fscanf(ifp, "%f", &Matrix_A[(i * order) + j]);

    /* Recebe os valores do arquivo para a matriz B */
    for(i = 0; i < order; i++)
        fscanf(ifp, "%f", &Vector_b[i]);
    
     
    float numerator; /* variavel utilizada para o calculo do erro */
    float denum;  /* variavel utilizada para o calculo do erro */
    numerator = 1;
    denum = 1;
    float x1;  /* variavel utilizada para o x(k+1) */
    int iter;  /* variavel utilizada para o calculo de iteracoes*/
    float maxl; /* variavel utilizada para o calculo da convergencia */
    float maxt;  /* variavel utilizada para o calculo do erro */
    int aux_iter=0; /* variavel auxiliar utilizada para o calculo de iteracoes*/
    
    /* Inicio do calculo da convergencia*/
    for(i = 0; i < order; i++)
    {
        maxl = 0;
        for(j = 0; j < order; j++)
        {
            if(i != j)
            {
              Matrix_A_e[(i * order) + j] =
                    Matrix_A[(i * order) + j] /
                    Matrix_A[(i * order) + i];
                maxl += Matrix_A_e[(i * order) + j];
            }
            else Matrix_A_e[(i * order) + j] = 0;
            maxt = max(maxt, maxl);
        }

        Vector_b_e[i] = Vector_b[i] /
            Matrix_A[(i * order) + i];
            
        Vector_x[i]=Vector_b_e[i];
    }

    if(maxt > 1)
        printf("Convergence criterion not satisfied.\n");
    else
        printf("Convergence criterion satisfied.\n");
    
    /****************************************************/
    /*************** METODO JACOBI-RICHARDSON ***********/
    /****************************************************/
    for(iter = 0; iter < ite_max;++iter){
        numerator = 0;
        denum = 0;
        aux_iter++;
        for(i = 0; i < order; i++){
            x1 = Vector_b_e[i];
            for (j = 0; j < order; j++) {
                if (i != j)
                  x1 -= Matrix_A_e[(i * order) + j] * Vector_x[j];
            }

            if (numerator < abs(x1 - Vector_x[i]) )
                numerator = abs(x1 - Vector_x[i]);
            if (denum < abs(x1))
                denum = abs(x1);
                
            Vector_x[i] = x1;

            if((numerator/denum) < (error) && numerator > (error*error)){
              iter = ite_max;}
        }
    }
      /* Calcula o resultado */  
    float res=0;  
    for(i = 0; i < order; i++)
        res += Matrix_A[(row_test * order) + i] * Vector_x[i];
        
    
    /* Imprime na tela os valores do resultado e iteraçao */    
    printf("-----------------------------------------\n");
    printf("Ite: %d\n",aux_iter);
    printf("RowTest: %d => [%f] =? [%f]\n",
            row_test, res, Vector_b[row_test]);
    printf("-----------------------------------------\n");

    
    
    

  
  return 0;
}